// ==================================================================================
// Shared Genomics Project MPI Codebase
// Version 1.0 30/04/2010
//
// (c) 2010 University of Manchester all rights reserved
//
// This file is distributed under the GNU General Public License, Version 2.  
// Please see the file COPYING.txt for more details
// ==================================================================================

/*!
\file
\brief Implementation of Quantitative QT Association (BitStrings)
*/

#include <math.h>

#include "passoc_qt.h"
#include "pmaths.h"

void passoc_qt_write_file_header(FILE *out) {
	if (out != NULL) fprintf(out, "CHR,BP,SNP,NMISS,BETA,SE,R2,T,P,EMP1,EMP2,BONF,HOLM,FDR_BH\n"); 
}

double passoc_qt_assoc(BitString ones_and_twos, struct sample *samples, int nSamples, struct small_locus *snp, FILE *out, double emp1, double emp2, struct adjusted *adj) {
	// Best model variables.
	register int i;
	double g_mean = 0, g_var = 0;
	double qt_mean = 0, qt_var = 0;
	int nanal = 0; // No. samples included in the analysis
	double qt_g_covar = 0;
	double beta, vbeta, t;

	// Iterate over individuals
	for (i = 0; i < nSamples; i++) {
		unsigned int mask;
		long bitpos = (long) i * 2;
		BOOL s1 = FALSE, s2 = FALSE;

		if (samples[i].pperson->missing) continue;

		// Read S1 value
		mask = 0x1;
		mask = mask << (bitpos % BITBULKSIZE);
		mask = ones_and_twos[bitpos/BITBULKSIZE] & mask;
		if (mask !=0) s1 = TRUE;
		bitpos++;

		// Read S2 value
		mask = 0x1;
		mask = mask << (bitpos % BITBULKSIZE);
		mask = ones_and_twos[bitpos/BITBULKSIZE] & mask;
		if (mask !=0) s2 = TRUE;

		if (!( s1 && (!s2))) {    //   10 = missing 
			qt_mean += samples[i].pperson->phenotype;
			if (!s1) {
				if (!s2)      //   00 = hom(11)
					g_mean+=2;
				else          //   01 = het(12)
					g_mean++;
			}
			nanal++;
		}
	}

	qt_mean /= (double) nanal;
	g_mean /= (double) nanal;

	// Iterate over individuals again
	for (i = 0; i < nSamples; i++) {
		unsigned int mask;
		long bitpos = (long) i * 2;
		BOOL s1 = FALSE, s2 = FALSE;

		if (samples[i].pperson->missing) continue;

		// Read S1 value
		mask = 0x1;
		mask = mask << (bitpos % BITBULKSIZE);
		mask = ones_and_twos[bitpos/BITBULKSIZE] & mask;
		if (mask !=0) s1 = TRUE;
		bitpos++;

		// Read S2 value
		mask = 0x1;
		mask = mask << (bitpos % BITBULKSIZE);
		mask = ones_and_twos[bitpos/BITBULKSIZE] & mask;
		if (mask !=0) s2 = TRUE;

		if ( !((s1) && (!s2))) {
			double g = 0;

			qt_var += 
				(samples[i].pperson->phenotype-qt_mean) * 
				(samples[i].pperson->phenotype - qt_mean);

			if (!s1) {
				if (!s2) g=2;
				else g=1;
			}

			g_var += (g-g_mean) * (g-g_mean) ;
			qt_g_covar += 
				(samples[i].pperson->phenotype - qt_mean) * 
				(g - g_mean) ;
		}
	}

	// Summary statistics
	qt_var /= (double) nanal - 1;
	g_var /= (double) nanal - 1;
	qt_g_covar /= (double) nanal - 1;

	// Test statistics
	beta = qt_g_covar / g_var;
	vbeta = ( qt_var/g_var - (qt_g_covar*qt_g_covar)/(g_var*g_var) ) / (nanal-2);
	t = beta / sqrt(vbeta);

	// Write-out to file if the FILE handle set.
	if (out > NULL && adj != NULL) {
		double t_p = pT(t, nanal - 2);
		double r2 =   (qt_g_covar * qt_g_covar ) / ( qt_var * g_var ) ;

		fprintf(out, "%i,%i,%s,%i,", snp->chr, snp->bp, snp->name, nanal);

		if (!realnum(beta)) fprintf(out, "NA,NA,NA,");
		else fprintf(out, "%f,%f,%f,", beta, sqrt(vbeta), r2);

		if (t_p >= 0) fprintf(out, "%f,%f,",  t, t_p);
		else  fprintf(out, "NA,NA,");

		// EMP1
		fprintf(out, "%f,%f,%f,%f,%f\n", emp1, emp2, adj->bonf, adj->holm, adj->fdr_bh);
	}

	// Capture the CHISQ for the calculation to find FDR et cetera.
	if (adj != NULL && out == NULL) {
		if (realnum(t)) {
			adj->chisq = t * t;
			adj->tcnt = nanal - 2;
		}
	}

	return t * t; // return chi-sq
}